R Notebook

source('../settings/settings.R')
source('commonFunctions.R')

persons <- SELECTED_SUBJECTS

all_Drive3 <- read.csv('..//data/processed/distancewise/TT1_Drive_3_30m_30m.csv')
all_Drive3$Subject <- as.factor(all_Drive3$Subject)
all_Drive3$logPerspiration <- log(all_Drive3$Perspiration)

mean_pp <- vector(mode="list", length=length(persons)) 
names(mean_pp) <- persons

std_pp <- vector(mode="list", length=length(persons)) 
names(std_pp) <- persons

# Mean (Turning)
mean_pp_seg0 <- vector(mode="list", length=length(persons)) 
names(mean_pp_seg0) <- persons
mean_pp_seg0_1 <- vector(mode="list", length=length(persons)) 
names(mean_pp_seg0_1) <- persons
mean_pp_seg0_2 <- vector(mode="list", length=length(persons)) 
names(mean_pp_seg0_2) <- persons
mean_pp_seg0_3 <- vector(mode="list", length=length(persons)) 
names(mean_pp_seg0_3) <- persons
mean_pp_seg0_4 <- vector(mode="list", length=length(persons)) 
names(mean_pp_seg0_4) <- persons

# Mean (Straight)
mean_pp_seg1 <- vector(mode="list", length=length(persons)) 
names(mean_pp_seg1) <- persons
mean_pp_seg2 <- vector(mode="list", length=length(persons)) 
names(mean_pp_seg2) <- persons
mean_pp_seg3 <- vector(mode="list", length=length(persons)) 
names(mean_pp_seg3) <- persons
mean_pp_seg4 <- vector(mode="list", length=length(persons)) 
names(mean_pp_seg4) <- persons
mean_pp_max <- vector(mode="list", length=length(persons)) 
names(mean_pp_max) <- persons

# SD (Turning)
std_pp_seg0 <- vector(mode="list", length=length(persons)) 
names(std_pp_seg0) <- persons
std_pp_seg0_1 <- vector(mode="list", length=length(persons)) 
names(std_pp_seg0_1) <- persons
std_pp_seg0_2 <- vector(mode="list", length=length(persons)) 
names(std_pp_seg0_2) <- persons
std_pp_seg0_3 <- vector(mode="list", length=length(persons)) 
names(std_pp_seg0_3) <- persons
std_pp_seg0_4 <- vector(mode="list", length=length(persons)) 
names(std_pp_seg0_4) <- persons

# SD (Straight)
std_pp_seg1 <- vector(mode="list", length=length(persons)) 
names(std_pp_seg1) <- persons
std_pp_seg2 <- vector(mode="list", length=length(persons)) 
names(std_pp_seg2) <- persons
std_pp_seg3 <- vector(mode="list", length=length(persons)) 
names(std_pp_seg3) <- persons
std_pp_seg4 <- vector(mode="list", length=length(persons)) 
names(std_pp_seg4) <- persons
std_pp_max <- vector(mode="list", length=length(persons)) 
names(std_pp_max) <- persons

for(p in persons) {
  pData <- all_Drive3[(all_Drive3$Subject==as.integer(p) | all_Drive3$Subject==p),]
  pData_act3 <- pData[pData$Activity==3,]
  
  pData_seg0 <- pData[pData$Phase==0,]
  pData_seg0_1 <- pData[pData$Activity==0 & pData$Phase==0 & pData$Time > 100 & pData$Time < 200,]
  pData_seg0_2 <- pData[pData$Activity==0 & pData$Phase==0 & pData$Time > 200 & pData$Time < 300,]
  pData_seg0_3 <- pData[pData$Activity==0 & pData$Phase==0 & pData$Time > 300 & pData$Time < 400,]
  pData_seg0_4 <- pData[pData$Activity==0 & pData$Phase==0 & pData$Time > 400 & pData$Time < 500,]
  
  pData_seg1 <- pData[pData$Phase==1 & pData$Activity==3 & pData$Time < 110,]
  pData_seg2 <- pData[pData$Phase==2 & pData$Activity==3 & pData$Time < 250,]
  pData_seg3 <- pData[pData$Phase==3 & pData$Activity==3 & pData$Time < 350,]
  pData_seg4 <- pData[pData$Phase==4 & pData$Activity==3,]
  
  mean_pp[[p]] <- mean(pData_act3$ppLogNormalized)
  std_pp[[p]] <- sd(pData$ppLogNormalized)
  
  mean_pp_seg0[[p]] <- mean(pData_seg0$ppLogNormalized)
  mean_pp_seg0_1[[p]] <- mean(pData_seg0_1$ppLogNormalized)
  mean_pp_seg0_2[[p]] <- mean(pData_seg0_2$ppLogNormalized)
  mean_pp_seg0_3[[p]] <- mean(pData_seg0_3$ppLogNormalized)
  mean_pp_seg0_4[[p]] <- mean(pData_seg0_4$ppLogNormalized)
  
  mean_pp_seg1[[p]] <- mean(pData_seg1$ppLogNormalized)
  mean_pp_seg2[[p]] <- mean(pData_seg2$ppLogNormalized)
  mean_pp_seg3[[p]] <- mean(pData_seg3$ppLogNormalized)
  mean_pp_seg4[[p]] <- mean(pData_seg4$ppLogNormalized)
  mean_pp_max[[p]] <- max(mean_pp_seg1[[p]], mean_pp_seg2[[p]], mean_pp_seg3[[p]], mean_pp_seg4[[p]])
  
  std_pp_seg0[[p]] <- sd(pData_seg0$ppLogNormalized)
  std_pp_seg0_1[[p]] <- sd(pData_seg0_1$ppLogNormalized)
  std_pp_seg0_2[[p]] <- sd(pData_seg0_2$ppLogNormalized)
  std_pp_seg0_3[[p]] <- sd(pData_seg0_3$ppLogNormalized)
  std_pp_seg0_4[[p]] <- sd(pData_seg0_4$ppLogNormalized)
  
  std_pp_seg1[[p]] <- sd(pData_seg1$ppLogNormalized)
  std_pp_seg2[[p]] <- sd(pData_seg2$ppLogNormalized)
  std_pp_seg3[[p]] <- sd(pData_seg3$ppLogNormalized)
  std_pp_seg4[[p]] <- sd(pData_seg4$ppLogNormalized)
  std_pp_max[[p]] <- max(std_pp_seg1[[p]], std_pp_seg2[[p]], std_pp_seg3[[p]], std_pp_seg4[[p]])
}

plt_AllAcc <- vector(mode="list", length=length(persons)) 
names(plt_AllAcc) <- persons

COLOR_ACC = "#02A3C8"
COLOR_PP = "#F28E8E"
COLOR_BRAKE = "#888888"

y1 <- list(
  tickfont = list(color = COLOR_ACC),
  title="Degree",
  range=c(0, max(all_Drive3$Acceleration))
)
y2 <- list(
  tickfont = list(color = COLOR_PP),
  overlaying = "y",
  side = "right",
  title = "Log Perspiration",
  showgrid = FALSE,
  range=c(-0.6, 0.9)
  # range=c(min(all_Drive3$ppLogNormalized), max(all_Drive3$ppLogNormalized))
)

for (p in persons) {
  pData <- all_Drive3[(all_Drive3$Subject==as.integer(p) | all_Drive3$Subject==p),]
  pData_seg0 <- pData[pData$Phase==0,]
  pData_seg1 <- pData[pData$Phase==1 & pData$Activity==3 & pData$Time < 110,]
  pData_seg2 <- pData[pData$Phase==2 & pData$Activity==3 & pData$Time < 250,]
  pData_seg3 <- pData[pData$Phase==3 & pData$Activity==3 & pData$Time < 350,]
  pData_seg4 <- pData[pData$Phase==4 & pData$Activity==3,]
  
  plot_Acc <- plot_ly(pData, x = ~Time, height=400, width=900) %>%
    # add_trace(name="Acceleration", y = ~Acceleration, type = 'scatter', mode = 'lines', line=list(width=1.5, color=COLOR_ACC)) %>% 
    add_trace(name="PP", y = ~ppLogNormalized, type = 'scatter', mode = 'lines', line=list(width=1.5, color=COLOR_PP), yaxis = "y2") %>%
    add_segments(x = min(pData$Time), xend = max(pData$Time), y = mean_pp[[p]], yend = mean_pp[[p]], 
                           yaxis = "y2", name="Avg. PP (straight)",
                           line=list(color="darkgray", dash = 'dot')) %>%
    add_segments(x = min(pData$Time), xend = max(pData$Time), y = mean_pp_seg0[[p]], yend = mean_pp_seg0[[p]], 
                           yaxis = "y2", name="Avg. PP (turning)",
                           line=list(color="black", dash = 'dot')) %>%
    add_segments(x = min(pData_seg1$Time), xend = max(pData_seg1$Time), y = mean_pp_seg1[[p]], yend = mean_pp_seg1[[p]], 
                           yaxis = "y2", name="Avg. PP (1st part)",
                           line=list(color="red", dash = 'dot')) %>%
    add_segments(x = min(pData_seg2$Time), xend = max(pData_seg2$Time), y = mean_pp_seg2[[p]], yend = mean_pp_seg2[[p]], 
                           yaxis = "y2", name="Avg. PP (2nd part)",
                           line=list(color="green", dash = 'dot')) %>%
    add_segments(x = min(pData_seg3$Time), xend = max(pData_seg3$Time), y = mean_pp_seg3[[p]], yend = mean_pp_seg3[[p]], 
                           yaxis = "y2", name="Avg. PP (3rd part)",
                           line=list(color="blue", dash = 'dot')) %>%
    add_segments(x = min(pData_seg4$Time), xend = max(pData_seg4$Time), y = mean_pp_seg4[[p]], yend = mean_pp_seg4[[p]], 
                           yaxis = "y2", name="Avg. PP (4th part)",
                           line=list(color="purple", dash = 'dot')) %>%
    layout(
      title=paste0("Subject #", p), 
      xaxis=list(title="Time [s]", range=c(0)), 
      yaxis=y1, 
      yaxis2=y2, 
      margin = list(l = 50, r = 50, b = 50, t = 50, pad = 4),
      legend = list(x = 0.5, xanchor = "center", y = 0.2, bgcolor = "rgba(0,0,0,0)", title="Metric", orientation = "h"),
      autosize = F
    )
  
  plt_AllAcc[[p]] <- plot_Acc
}


htmltools::tagList(plt_AllAcc)

NUMBER_OF_CLUSTERS = 3

color_darkpink = "#e75480"
CLUSTER_BRANCH_COLORS <- c("blue", "darkred", color_darkpink, "black")[1:NUMBER_OF_CLUSTERS]
CLUSTER_LABEL_COLORS <- c("blue", "darkred", color_darkpink, "black")[1:NUMBER_OF_CLUSTERS]

dfPP <- as.data.frame(cbind(
                            unlist(mean_pp), 
                            unlist(std_pp), 
                            unlist(mean_pp_seg0), 
                            unlist(mean_pp_seg1), 
                            unlist(mean_pp_seg2), 
                            unlist(mean_pp_seg3), 
                            unlist(mean_pp_seg4),
                            unlist(mean_pp_max),
                            unlist(std_pp_seg0), 
                            unlist(std_pp_seg1), 
                            unlist(std_pp_seg2), 
                            unlist(std_pp_seg3), 
                            unlist(std_pp_seg4),
                            unlist(std_pp_max),
                            unlist(mean_pp_seg0_1),
                            unlist(mean_pp_seg0_2),
                            unlist(mean_pp_seg0_3),
                            unlist(mean_pp_seg0_4),
                            unlist(std_pp_seg0_1),
                            unlist(std_pp_seg0_2),
                            unlist(std_pp_seg0_3),
                            unlist(std_pp_seg0_4)
                    ))

names(dfPP) <- c("MeanPP", "StdPP", 
                 "MeanPP_Seg0", "MeanPP_Seg1", "MeanPP_Seg2", "MeanPP_Seg3", "MeanPP_Seg4", "MeanPP_SegMax",
                 "StdPP_Seg0", "StdPP_Seg1", "StdPP_Seg2", "StdPP_Seg3", "StdPP_Seg4", "StdPP_SegMax",
                 "MeanPP_Seg0_1", "MeanPP_Seg0_2", "MeanPP_Seg0_3", "MeanPP_Seg0_4",
                 "StdPP_Seg0_1", "StdPP_Seg0_2", "StdPP_Seg0_3", "StdPP_Seg0_4")

behavioralMatrixClustering <- as.matrix(dfPP)

distMatrix <- dist(behavioralMatrixClustering)
hresults <- distMatrix %>% hclust(method="average")

hc <- hresults %>% 
      as.dendrogram %>%
      set("nodes_cex", NUMBER_OF_CLUSTERS) %>%
      set("labels_col", value = CLUSTER_LABEL_COLORS, k=NUMBER_OF_CLUSTERS) %>%
      # set("leaves_pch", 19) %>%
      # set("leaves_col", value = c("gray"), k=NUMBER_OF_CLUSTERS) %>%    
      set("branches_k_color", value=CLUSTER_BRANCH_COLORS, k=NUMBER_OF_CLUSTERS)

plot(hc)
legend("topright", 
     title="Drive=Motoric \nHierachical Clustering",
     legend = c("Group 1", "Group 2", "Group 3"), 
     col = c("darkred", "pink" , "blue"),
     pch = c(20,20,20), bty = "n",  pt.cex = 1.5, cex = 0.8 , 
     text.col = "black", horiz = FALSE, inset = c(0.4, 0.1))

Clustering with SD of PP

library(dendextend)

clusteringDf <- dfPP %>% select(StdPP)
NUMBER_OF_CLUSTERS = 3

color_darkpink = "#e75480"
CLUSTER_BRANCH_COLORS <- c("darkred", "blue", "red")[1:NUMBER_OF_CLUSTERS]
CLUSTER_LABEL_COLORS <- c("darkred", "blue", "red")[1:NUMBER_OF_CLUSTERS]

behavioralMatrixClustering <- as.matrix(clusteringDf)
rownames(behavioralMatrixClustering) <- paste0("#", persons)
distMatrix <- dist(behavioralMatrixClustering)
hresults <- distMatrix %>% hclust

hc <- hresults %>% 
      as.dendrogram %>%
      set("nodes_cex", NUMBER_OF_CLUSTERS) %>%
      set("labels_col", value = CLUSTER_LABEL_COLORS, k=NUMBER_OF_CLUSTERS) %>%
      # set("leaves_pch", 19) %>%
      # set("leaves_col", value = c("gray"), k=NUMBER_OF_CLUSTERS) %>%    
      set("branches_k_color", value=CLUSTER_BRANCH_COLORS, k=NUMBER_OF_CLUSTERS)

plot(hc)
legend("topright", 
     title="Drive=Motoric \nChange of Arousal",
     legend = c("Exceptional SD" , "Low SD" , "High SD"), 
     col = c("darkred", "blue" , "red"),
     pch = c(20,20,20), bty = "n",  pt.cex = 1.5, cex = 0.8 , 
     text.col = "black", horiz = FALSE, inset = c(0.0, 0.1))

library(cluster)
silhouette_score <- function(k){
  km <- kmeans(clusteringDf, centers = k, nstart=25)
  ss <- silhouette(km$cluster, dist(clusteringDf))
  mean(ss[, 3])
}
k <- 2:10
avg_sil <- sapply(k, silhouette_score)
plot(k, type='b', avg_sil, xlab='Number of clusters', ylab='Average Silhouette Scores', frame=FALSE)

# Store clustering data
fPath <- str_interp("../data/processed/analysis/TT1_Drive_3_PP.csv")
dfx <- dfPP
dfx <- cbind(persons, dfx)
names(dfx) <- c("Subject" , names(dfPP))
write.csv(dfx, fPath, row.names = F)